
\chapter{常微分方程}
常微分方程，即一元函数的微分方程；偏微分方程，就是多元函数的微分方程。
在这一章，我们先介绍初值问题，然后介绍边值问题。
\section{初值问题}
如果函数$y$在定义域一端$x=a$处的信息已知，需要推测其他地方的函数值，那么这种问题就叫做初值问题。
一个典型的例子是，将单摆拉开一定角度，然后松手，推测一天以后单摆在哪里，速度是多少。
将松手的瞬间定位$t=0$，则
\begin{eqnarray}
\theta(t=0) &=& \theta_0, \\
\dot{\theta}(t=0) &=& 0
\end{eqnarray}
定义了初始时刻的摆角和角速度。
然后根据牛顿第二定律做动力学的分析，推测$\theta(t)$与$\dot{\theta}(t)$的轨迹，这就是初值问题。

相对的，边值问题是这样的一类问题，我们知道定义域一端$x=a$处$y$的一部分信息，以及另一端$x=b$的一部分信息，然后根据动力学方程推断出整个定义域内各处的$y$值。
可怜的中学生们也可能遇到这样的问题，比如，$t=0$时刻，单摆的摆角为
\begin{equation}
\theta(t=0) = \theta_1,
\end{equation}
但是不知道它的角速度是多少。
而在$t=100s$，单摆的摆角为
\begin{equation}
\theta(t=10) = \theta_2,
\end{equation}
也不知道角速度是多少，求单摆的最大摆角$\theta_0$。
这就是一个时间轴上的边值问题，在下一章我们介绍相关的解法。

\subsection{微分方程降阶}
如果方程中最高阶的导数是$n$阶导数$y^{(n)}$，我们就把这个方程叫做$n$阶微分方程，
\begin{equation}
\frac{d^n}{dx^n} y = f(y, \frac{dy}{dx}, \cdots, \frac{d^{n-1}}{dx^{n-1}}y),
\end{equation}
一般需要$n-1$个关于$y,y',y'',\cdots$的初始条件。
为了方便处理，我们一般把高阶的微分方程改写为一阶微分方程组，定义$n$个函数$y^{0}, \cdots, y^{(n-1)}$，分别表示函数$y$的$0,\cdots,n-1$阶导函数，则原来的微分方程可以写作，
\begin{eqnarray}
\frac{d}{dx} y^{(0)} &=& y^{(1)}, \\
\frac{d}{dx} y^{(1)} &=& y^{(2)}, \\
& \vdots &  \\
\frac{d}{dx} y^{(n-1)} &=& f(y^0, \cdots, y^{(n-1)}), \\
\end{eqnarray}
这$n$个1阶微分方程，加上$n-1$个关于$y^{(0)},y^{(1)},\cdots,y^{(n-1)}$的初始条件，唯一地确定了解。

为了简单形象，我们举一个例子，一维谐振子，即光滑平面上固定在横向理想弹簧上的滑块，本来在平衡位置，现在拉开一段距离$a$（弹性形变），然后松开，谐振子的运动方程为
\begin{equation}
m \frac{d^2}{dt^2} x(t) = - kx(t),
\end{equation}
为了简便，我们设$m=k=a=1$（也可以通过定义新的变量将方程变为无量纲方程，效果是一样，简便）。
那么方程变为
\begin{equation}
\frac{d^2}{dt^2} x = -x,
\end{equation}
初始条件为$x(0)=1, x'(0)=0$，为了降阶，我们定义$v(t)=\frac{dx}{dt}$，方程变为
\begin{eqnarray}
\frac{dx}{dt} &=& v, \\
\frac{dv}{dt} &=& - x.
\end{eqnarray}


\subsection{欧拉法}
我们的目标是求得$y(x)$在定义域中不同地方的值。
于是我们将定义域$\left[x_0=a, x_n=b\right]$等分为$n$份，步长为$h=(b-a)/n$，得到离散点$x_0, x_1, x_2, \cdots, x_{n-1}, x_n$，然后求取$f(x_i)$的近似值$y_i = f(x_i)$，就算是完成了数值任务——想象计算机高兴地说“我解出来了！”。

最简单直接的是欧拉法，做$y(x)$的一阶泰勒展开，
\begin{eqnarray}
y(x_{i+1}) = y(x_i) + h \frac{dy}{dx}|_{x_i} + O(h^2) = y(x_i)+ hf(x_i, y(x_i)) + O(h^2),
\end{eqnarray}
所以可以构造迭代公式
\begin{eqnarray}
y_{i+1} = y_i + hf(x_i, y_i)
\end{eqnarray}
完成 $y_i$ 到 $y_{i+1}$ 的迭代。
注意，上面我们用$y(x_i)$代表函数值的理论值，用$y_i$代表$y(x_i)$的数值近似。
在每一步迭代，截断误差约为 $O(h^2)$，所以迭代求得 $y_n$ 时，总误差 $y_n - y(x_n) \propto nh^2 \propto h$，所以欧拉法的精度为$O(h)$，即一阶精度。

我们可以想到一个很简单的例子，即雨滴下落。
雨滴从云中掉落，重力使它持续加速，如果没有空气阻力，就会像水做的子弹一样，把人类砸伤。
还好有空气阻力，空气阻力与雨滴速度的平方成正比，使得雨滴速度增加到某个值，就不再增加了。
所以雨滴落地时速度一般是几米每秒，人类可以愉快地在雨中歌唱，并写下著名舞曲“Singing in the rain.”。

在这个故事中，雨滴速度的初始值是
\begin{equation}
v(0) = 0,
\end{equation}
其微分方程为
\begin{equation}
\frac{dv}{dt} = g - \alpha v^2,
\end{equation}
其中$\alpha$是参数，与雨滴质量等因素有关，我们简单将它设为1，也许不准确，但不妨碍我们展示这个小故事。
我们用欧拉法解微分方程，确定$v(t)$，代码如下
\lstset{basicstyle=\footnotesize}
\begin{lstlisting}

#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>

/* Euler(...)           solve 1-order ODE for initial problem
 * double a             [a,b]
 * double b
 * double *func(double x, double y)     f(x,y)=dy/dx
 * int n                [a,b] is divided into n pieces
 * double y0            initial value y(a)
 * double *y            y[n+1] stores y0, y1,..., yn
 */
void Euler(double a, double b, double (*func)(double x, double y), int n, double y0, double *y){
        
        int i;
        double h = (b-a)/n; //步长
        y[0] = y0;
        for(i=0;i<n;i++){
                y[i+1] = y[i] + h*func(a + i*h, y[i]);
        }
        ofstream fp("Euler.txt"); //将结果输出到Euler.txt，用作画图
        for(i=0;i<=n;i++){
                fp<<a+i*h<<"\t"<<y[i]<<endl;
        }
        fp.close();
}

double f(double x, double y){ // dv/dt = 10 - v^2
        return 10 - y*y;
}

int main(){
        
        int i,j,k;
        
        int n=10000;
        double *y = new double [n];
        Euler(0, 2, f, n, 0, y);
        
        delete [] y;
        return 0;
}       
\end{lstlisting}
跑完程序以后，我们得到程序文件夹下的Euler.txt文件，就可以用gnuplot画图，gnuplot定位到当前文件夹后，打入
\begin{lstlisting}
plot "Euler.txt" w l
\end{lstlisting}
即用折线连接Euler.txt文件中给出的离散点，绘成图片(\ref{fig:Euler-raindrop}),图中我们可以看到，雨滴速度变大，然后趋于饱和的过程。
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{Euler}
\caption{雨滴速度(m/s) v.s. 时间(s)。}
\label{fig:Euler-raindrop}
\end{figure}

\subsection{改进的欧拉法}
对上面描述的这种初值问题，$y(x_{i+1})$的积分形式为
\begin{equation}
y(x_{i+1}) = y(x_i) + \int_{x_i}^{x_{i+1}} f(x,y) dx,
\end{equation}
上式右边的积分，如果做如下矩形近似
\begin{eqnarray}
\int_{x_i}^{x_{i+1}} f(x,y) dx &=& \int_{x_i}^{x_{i+1}} (f(x_i,y(x_i)) + O(x) ) dx   \nonumber \\
&\approx& h f(x_i,y(x_i)) + O(h^2),
\end{eqnarray}
则得到欧拉法。

那么如果我们希望得到更好的近似，可以用梯形法构造上述积分，
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = \frac{h}{2}[ f(x_i,y(x_i)) + f(x_{i+1}, y(x_{i+1}) ) ] + O(h^3),
\end{equation}
但是我们不知道$y(x_{i+1})$，所以做近似
\begin{equation}
y(x_{i+1}) = y(x_i) + h f(x_i,y_i) + O(h^2),
\end{equation}
得到
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = \frac{h}{2}[ f(x_i,y(x_i)) + f(x_{i+1}, y(x_i) + h f(x_i,y_i)) ]  + O(h^3),
\end{equation}
用相应的数值近似代替上式中的解析表达式，得到迭代公式
\begin{equation}
y_{i+1} = y_i + \frac{h}{2}[ f(x_i,y_i) + f(x_{i+1}, y_i + h f(x_i,y_i)) ],
\end{equation}
这就是改进的欧拉法，每步具有$O(h^3)$精度，推到$y_n$时，具有$O(h^2)$精度。

\subsection{龙格库塔方法}
延续上面的思路，如果用中值法处理积分，则得到二阶龙格库塔方法，
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = h f(x_i + h/2, y(x_i+h/2)) + O(h^3),
\end{equation}
未知的$f(x_i + h/2, y(x_i+h/2))$可以用下面的近似公式替代，
\begin{eqnarray}
f(x_i + \frac{h}{2}, y(x_i+h/2)) &=& f(x_i+\frac{h}{2}, y(x_i) + \frac{h}{2}f(x_i,y(x_i)) + O(h^2)) \nonumber\\
&=& f(x_i+\frac{h}{2}, y(x_i) + \frac{h}{2}f(x_i,y(x_i))) + O(h^2),
\end{eqnarray}
得到公式
\begin{equation}
\int_{x_i}^{x_{i+1}} f(x,y) dx = h f(x_i + \frac{h}{2}, y(x_i) + \frac{h}{2}f(x_i, y(x_i)) ) + O(h^3),
\end{equation}
相应的，数值迭代公式为
\begin{equation}
y_{i+1} = y_i + h f(x_i + \frac{h}{2}, y_i + \frac{h}{2} f(x_i,y_i)),
\end{equation}
或者写得稍微清晰一些
\begin{equation}
\left\{
\begin{aligned}
&K_1 = f(x_i, y_i)\\
&y_{i+1} = y_i + h f(x_i + h/2, y_i + (h/2)K_1)
\end{aligned}
\right.
\end{equation}
这个公式每步具有$O(h^3)$，推到$y_n$时累积有$O(h^2)$精度，所以叫做二阶龙格库塔公式。

若需要四阶龙格库塔公式，即需要每步精度为$O(h^5)$，我们想起了辛普生方法，即
\begin{equation}
\int^{x_i +h}_{x_i} f(x,y) dx = \frac{h}{6}[ f(x_i,y(x_i)) + 4f(x_i + \frac{h}{2},y(x_i+\frac{h}{2})) + f(x_i + h,y(x_i+h)) ] + O(h^5),
\end{equation}
找到$f(x_i + \frac{h}{2},y(x_i+\frac{h}{2}))$，$f(x_i + h,y(x_i+h))$的4阶精度替代方程，即可构造出4阶龙格库塔公式，下面的数值迭代式是前人的构造，
\begin{eqnarray}
\left\{
\begin{aligned}
&K_1 = f(x_i, y_i)\\
&K_2 = f(x_i + h/2, y_i + (h/2)K_1)\\
&K_3 = f(x_i + h/2, y_i + (h/2)K_2)\\
&K_4 = f(x_i + h, y_i + hK_3)\\
&y_{i+1} = y_i + \frac{h}{6}(K_1 + 2K_2 + 2K_3 + K_4)
\end{aligned}
\right.
\end{eqnarray}
这个公式在每个步长的误差为$O(h^5)$，累积误差为$O(h^4)$。
上面这个式子中，$K1$是$x_i$处的斜率，而$K2,K3$是$x_i+h/2$处斜率的两个近似，$K4$是$x_i+h$处的斜率。
所以四阶龙格库塔方法利用了3处的斜率信息，然后进行加权，所以能得到对$y_{i+1}$的更精确的推断。
四阶龙格库塔的应用相当普遍，下面的代码实现一个步长的计算
\begin{lstlisting}
/* rk4:                 4阶龙格库塔一个步长的计算
 * double x:            xi
 * double h:            步长
 * double *y:           xi处已知的 y, y', y'', ..., y^{(n-1)}
 * double *ynext        计算的结果 xi+h 处的 y, y', y'', ..., y^{(n-1)}
 * void (*deriv)(int n, double x, double *y, double *dydx)
 *                      xi 处的导函数，根据微分方程分别计算 y, y', ..., y^{(n-1)} 在 x 处的导函数，并储存在 double *dydx 中
 */
void rk4(double x, double h, int n, double *y, double *ynext, void (*deriv)(int n, double x, double *y, double *dydx)){
        int i,j;
        double hh = h*0.5;
        double h6 = h/6;
        double *K, *ytemp;
        K = new double [n];
        deriv(n, x, y, K);// K1
        for(i=0;i<n;i++) ynext[i] = K[i];

        ytemp = new double [n];
        for(i=0;i<n;i++) ytemp[i] = y[i] + hh*K[i];// yi + h/2*K1
        deriv(n, x+hh, ytemp, K);// K2 
        for(i=0;i<n;i++) ynext[i] += 2 * K[i];

        for(i=0;i<n;i++) ytemp[i] = y[i] + hh*K[i];// yi + h/2*K2
        deriv(n, x+hh, ytemp, K);// K3 
        for(i=0;i<n;i++) ynext[i] += 2 * K[i];

        for(i=0;i<n;i++) ytemp[i] = y[i] + h*K[i];// yi + h*K3
        deriv(n, x+h, ytemp, K);// K4
        for(i=0;i<n;i++) ynext[i] += K[i];

        for(i=0;i<n;i++){// update y_{j+1}
                ynext[i] = y[i] + h6*ynext[i];
        }
        delete [] K;
        delete [] ytemp;
}
\end{lstlisting}

\subsection{经典谐振子}
下面演示怎样用4阶龙格库塔方法，求解经典谐振子，前面已经提到，谐振子的运动方程为
\begin{equation}
m \frac{d^2}{dt^2} x(t) = - kx(t),
\end{equation}
为了简便，我们设$m=k=a=1$，则方程变为
\begin{equation}
\frac{d^2}{dt^2} x = -x,
\end{equation}
初始条件为$x(0)=1, x'(0)=0$，则有
\begin{eqnarray}
\frac{dx}{dt} &=& v, \\
\frac{dv}{dt} &=& - x.
\end{eqnarray}
然后利用龙格库塔方法，即可实现求解，代码如下
\begin{lstlisting}
头文件声明。。。

void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];// dx/dt = v
        dydx[1] = - y[0];// dv/dt = -x
}

int main(){
        int i,j, n=2, N=1000;
        double a=0, b=10;
        double h = (b-a)/N;
        double *y = new double [n];
        double *ynext = new double [n];
        ofstream fp("classic_HO.txt");//输出文件流对象，将每步求得的结果输出到文件
        y[0]=1; y[1] = 0;
        fp<<a<<"\t"<<y[0]<<"\t"<<y[1]<<"\n";
        for(i=0;i<N;i++){
                rk4(a+i*h, h, n, y, ynext, derivative);
                for(j=0;j<n;j++) y[j] = ynext[j];
                fp<<a+(i+1)*h<<"\t"<<ynext[0]<<"\t"<<ynext[1]<<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
得到的数据存放在"classic\_Ho.txt"中，然后可以用gnuplot画图
\begin{lstlisting}
plot "classic_HO.txt" u 1:2 w l lw 3 t "x (m), Arial, 15", "" u 1:3 w l lw 3 t "v (m/s)"
\end{lstlisting}
上面的gnuplot指令表示，利用文件中的第1、2列分别作为横、纵坐标画一条曲线，第1、3列分别作为横、纵坐标画另一条曲线。
得到的图片如图\ref{fig:ClassicHO},
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{Classic_HO}
\label{fig:ClassicHO}
\caption{经典谐振子：4阶龙格库塔求解。}
\end{figure}

\paragraph{阻尼振动}
如果谐振子浸泡在液体中，它的运动就会受到液体的阻力，运动速度较小时，阻力与它的速度成正比，方向相反，
\begin{equation}
F_f = - \gamma \frac{dx}{dt} = - \gamma v,
\end{equation}
所以振子的运动方程变为
\begin{equation}
m \frac{d^2}{dt^2} x = - k x - \gamma v,
\end{equation}
因为我们假设$m = k = 1$，所以方程变为
\begin{equation}
\frac{d^2}{dt^2} x = - x - \gamma v,
\end{equation}
方程降阶以后得到
\begin{eqnarray}
\frac{d}{dt} x &=& v, \\
\frac{d}{dt} v &=& -x - \gamma v.
\end{eqnarray}
我们只需要稍微修改上面代码的$derivative$函数，即可计算阻尼运动。
\begin{lstlisting}
void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];// dx/dt = v
        dydx[1] = - y[0] - 0.1 * y[1];// dv/dt = -x - gamma v
}
\end{lstlisting}
为了看得更清楚，我们计算$t=[0,100]s$的振动，得到图\ref{fig:HOdamping}中所示的结果。
可以看到，因为液体阻尼的存在，简谐振动的振幅快速衰减。
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{classicHO_damping}
\label{fig:HOdamping}
\caption{谐振子阻尼运动，$m=k=a=1$，$\gamma=0.1$。}
\end{figure}

\paragraph{共振}
如果除了阻尼以外，谐振子还受一个周期性的驱动力，
\begin{equation}
F = F_0 \sin(\omega t),
\end{equation}
那么谐振子的运动方程就变为
\begin{equation}
m \frac{d^2}{dt^2} x = - k x - \gamma v + F_0 \sin(\omega t),
\end{equation}
如果我们假设$m=k=a=1$，则有
\begin{eqnarray}
\frac{d}{dt} x &=& v, \\
\frac{d}{dt} v &=& -x - \gamma v + F_0 \sin(\omega t).
\end{eqnarray}
我们设$\gamma=0.1, F_0 = 0.05$，分别取$\omega=0.8,0.9,1.0,1.1$，得到图\ref{fig:classicHO_Resonance}。
可以看到，外部驱动力的频率与本征频率相同，即$\omega=1.0$的时候，谐振子振幅最大。
\begin{figure}
\centering
\includegraphics[width=1.0\textwidth]{classicHO_Resonance}
\label{fig:classicHO_Resonance}
\caption{谐振子在外部驱动力下的共振现象。}
\end{figure}

\subsection{Numerov方法}
对于如下形式的常微分方程
\begin{eqnarray}
 \frac{ d^2 y}{dx^2} = - g(x) y(x) + s(x),
 \label{Numerov1}
\end{eqnarray}
其中 $g(x), s(x)$ 的形式是已知的，前苏联天文学家 Numerov 提出了5阶精度的算法：
\begin{eqnarray}
y(x_{n+1}) [ 1 + \frac{h^2}{12} g(x_{n+1}) ] =&& 2 y(x_n) [ 1 - \frac{ 5 h^2}{12} g(x_n) ] - y(x_{n-1}) [1 + \frac{h^2}{12} g(x_{n-1}) ]
\nonumber\\
&& + \frac{ h^2 }{12} [ s(x_{n+1}) + 10 s(x_n) + s(x_{n-1}) ] + O(h^6),
\end{eqnarray}
如果用节点上的计算值 $y_{n+1}, y_n, y_{n-1}$ 代替上式中的 $y(x_{n+1}), y(x_n), y(x_{n-1})$，并忽略上式右边最后一项 $O(h^6)$，则得到一个精度为$O(h^6)$的数值递推式
\begin{eqnarray}
y_{n+1} ( 1 + \frac{h^2}{12} g_{n+1}) &=& 2 y_n ( 1 - \frac{ 5 h^2}{12} g_n ) - y_{n-1} (1 + \frac{h^2}{12} g_{n-1} ) \nonumber\\
&& + \frac{ h^2 }{12} ( s_{n+1} + 10 s_n + s_{n-1} ),
\end{eqnarray}
其中 $g_n = g(x_n), ~ s_n = s(x_n)$。
如果区间最左边的两个节点函数值 $y_0, y_1$ 已知，就可以推得 $y_2, y_3, \cdots, y_n$，即得到方程的数值解。
\paragraph{推导}
下面我们推导 Numerov 方法的公式。
$y(x)$ 的泰勒展开式有
\begin{eqnarray}
y( x_{n+1} ) &=& y(x_n) + h y'(x_n) + \frac{ h^2}{2} y''(x_n) + \frac{ h^3}{6} y'''(x_n) 
\nonumber\\
&& + \frac{ h^4}{24} y''''(x_n) + \frac{ h^5}{120} y'''''(x_n) + O( h^6),
\nonumber\\
y( x_{n-1} ) &=& y(x_n) - h y'(x_n) + \frac{ h^2}{2} y''(x_n) - \frac{ h^3}{6} y'''(x_n) 
\nonumber\\
&& + \frac{ h^4}{24} y''''(x_n) - \frac{ h^5}{120} y'''''(x_n) + O( h^6),
\nonumber\\
\end{eqnarray}
上两式相加，得到
\begin{eqnarray}
y( x_{n+1} ) + y( x_{n-1} ) - 2 y(x_n) = h^2 y'' (x_n) + \frac{ h^4 }{12} y''''(x_n) + O(h^6).
\label{Numerov2}
\end{eqnarray}
结合（\ref{Numerov1}），可以推得
\begin{eqnarray}
y( x_{n+1} ) + y( x_{n-1} ) - 2 y(x_n) = h^2 ( - g(x_n) y(x_n) + s(x_n) ) + \frac{ h^4 }{12} y''''(x_n) + O(h^6).
\label{Numerov3}
\end{eqnarray}
$ y''''(x_n) $ 可以用（\ref{Numerov2}）类似的方法得到（用 $y''$ 代替式中的 $y$ 即得下式），
\begin{eqnarray}
y''(x_{n+1}) + y''(x_{n-1}) - 2 y''(x_n) = h^2 y'''' (x_n) + O(h^4), 
\end{eqnarray}
将（\ref{Numerov1}）代入上式得到，
\begin{eqnarray}
&& - g( x_{n+1} ) y( x_{n+1} ) + s( x_{n+1} ) - g( x_{n-1} ) y( x_{n-1} ) + s( x_{n-1} ) 
\nonumber\\
&& - 2 ( -g(x_n) y(x_n) + s(x_n) ) = h^2 y''''(x_n) + O(h^4), 
\end{eqnarray}
将上式代入（\ref{Numerov3}），并整理，即得
\begin{eqnarray}
y(x_{n+1}) [ 1 + \frac{h^2}{12} g(x_{n+1}) ] =&& 2 y(x_n) [ 1 - \frac{ 5 h^2}{12} g(x_n) ] - y(x_{n-1}) [1 + \frac{h^2}{12} g(x_{n-1}) ]
\nonumber\\
&& + \frac{ h^2 }{12} [ s(x_{n+1}) + 10 s(x_n) + s(x_{n-1}) ] + O(h^6).
\end{eqnarray}

\section{边值问题}
最简单的是两点边值问题，在一个边值处的条件不足以限制唯一解，需要结合另一边的条件才行。
为了直观形象，我们可以再次访问一维无限深势阱薛定谔方程。
\begin{equation}
V(x) = \left\{
\begin{aligned}
&0, ~~~ 0<x<a\\
&\infty, ~~else
\end{aligned}
\right.
\end{equation}
在$[0,a]$区间内，薛定谔方程为
\begin{equation}
- \frac{ \hbar^2 }{2m} \frac{ d^2 }{dx^2} \psi = E \psi，
\end{equation}
相应的本征值和解为
\begin{eqnarray}
E_n &=& \frac{ n^2 \pi^2 \hbar^2 }{ 2 m a^2}, \\
\psi_n &=& \sqrt{\frac{2}{a}} \sin (n \pi x /a  ).
\end{eqnarray}

在数值中，为了简便，我们取$\frac{\hbar^2}{2m} = a = 1$，则薛定谔方程变为
\begin{equation}
- \frac{d^2}{dx^2} \psi = E \psi, 
\end{equation}
其解析解为
\begin{eqnarray}
E_n &=& n^2 \pi^2 , \\
\psi_n &=& \sqrt{2} \sin (n \pi x  ).
\end{eqnarray}

在矩阵本征值一章，我们已经用有限差分将微分方程转化为矩阵问题，求解了这个问题。
在这里，我们将它作为边值问题的典型，用微分方程的方法求解。
在求解时，本征值是未知的，所以不能像初值问题中一样，从左推到右。

\subsection{打靶法}
我们可以换一个角度看问题，如果已知$\psi(0)=0, \psi'(0)=1$\footnote{这里$\psi'(0)$可以设为不同的值，仅与波函数归一化有关，不影响解的唯一性}，如果已知本征值$E$，则可以用龙格库塔方法从左到右推得$\psi_1$的值。
所以，换一句话说，这样做的话，$\psi_1$是本征值$E$的函数，我们要求的，是使得$\psi_1$为等于$0$的$E$值
\begin{equation}
\psi_1 (E) = 0,
\end{equation}
这样，问题转化为了方程求根问题，每个根对应一个本征值，有了具体的本征值，就可以用龙格库塔方法从左往右推出本征波函数。

如果用牛顿法求解上面的方程，就对应人们比较常说的打靶法。
如果对效率要求不高，也可以用更稳定的二分法。

要做龙格库塔，先把薛定谔方程降阶处理
\begin{eqnarray}
\frac{d}{dx} y^{(0)} &=& y^{(1)}, \\
\frac{d}{dx} y^{(1)} &=& - E y^{(0)}
\end{eqnarray}
然后在程序中定义$\psi_1 (E)$，然后做二分法即可，代码如下
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>
#include"library.h"

double E;//全局变量E

void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];
        dydx[1] = - E* y[0];
}

double psib(double x){

        E=x;//将x的值赋给全局变量E

        int N=10000,n=2;
        double a=0, ya=0, deriv_a=1, b=1;
        int i,j;
        double h = (b-a)/N;
        double *y = new double [n];
        double *ynext = new double [n];
        y[0]=ya; y[1] = deriv_a;
        for(i=0;i<N;i++){
                rk4(a+i*h, h, n, y, ynext, derivative);
                for(j=0;j<n;j++) y[j] = ynext[j];
        }
        return y[0];//返回 psi(b)
}

int main(){

        double eig = bisection(psib, 50, 100, 1E-9);
        cout<<"E/pi/pi="<<eig/M_PI/M_PI<<endl;

        return 0;
}
\end{lstlisting}

\subsection{格林函数法}
{\bf 格林定理}：如果$u,v$是标量函数，则有
\begin{eqnarray}
\vec{\nabla} \cdot ( u \vec{\nabla}v ) &=& u \nabla^2 v + (\nabla u) \cdot (\nabla v), \\
\vec{\nabla} \cdot ( v \vec{\nabla}u ) &=& v \nabla^2 u + (\nabla v) \cdot (\nabla u), 
\end{eqnarray}
所以有
\begin{equation}
\int_V (v \nabla^2 u - u \nabla^2 v) d \tau = \oint_S (v \vec{\nabla} u - u \vec{\nabla} v ) \cdot d \vec{\sigma}.
\end{equation}
其中$V$是一个封闭体积，$S$是其表面。
这个式子就是格林定理，我记得在电动力学唯一性定理的证明中就用到了。

{\bf 广义格林定理}：上面的定理很容易推广到如下情形，$L$是算子
\begin{equation}
L = \vec{\nabla} \cdot [ p(\vec{r}) \vec{\nabla}] + q(\vec{r}),
\end{equation}
则有
\begin{equation}
\int_V (v L u - u L v) d \tau = \oint_S p (v \vec{\nabla} u - u \vec{\nabla} v ) \cdot d \vec{\sigma}.
\end{equation}
格林定理是上式的$p(\vec{r})=1, q(\vec{r})=0$的特殊情形。

广义格林定理可以如下证明。
\begin{equation}
L u = \vec{\nabla} \cdot [ p(\vec{r}) \vec{\nabla} ] u + q(\vec{r}) u,
\end{equation}
将上式左乘$v$，然后对体积$V$做积分，得到
\begin{eqnarray}
\int_V vLu d\tau &=& \int_V ( v \vec{\nabla} \cdot [ p(\vec{r}) \vec{\nabla} ] u + q(\vec{r}) vu ) d \tau \\
&=& \int_V ( \vec{\nabla} \cdot (vp \vec{\nabla}u) - p \vec{\nabla}v \cdot \vec{\nabla}u + q(\vec{r}) vu ) d \tau,
\end{eqnarray}
类似地，可以写出
\begin{equation}
\int_V uLv d\tau = \int_V ( \vec{\nabla} \cdot (up \vec{\nabla}v) - p (\vec{\nabla}u) \cdot (\vec{\nabla}v) + q(\vec{r}) vu ) d \tau,
\end{equation}
两式相减，得到
\begin{equation}
\int_V (vLu - uLv) d \tau = \oint_S p(v \vec{\nabla}u - u \vec{\nabla} v) \cdot d\vec{\sigma}
\end{equation}
这就证明了广义格林定理。

很多物理问题中，当$r \rightarrow \infty$时，$pv \vec{\nabla}u, pu \vec{\nabla} v$都按$\frac{1}{r}$的2以上阶数衰减，所以上式右侧在无穷远处趋于0。
于是我们会得到
\begin{equation}
\int_V vLu d\tau = \int_V uLv d \tau.
\end{equation}
那么，如果待求的2阶偏微分方程为
\begin{equation}
L u = f,
\end{equation}
而我们设计出了$v(\vec{r}, \vec{r'})$，使得
\begin{equation}
L v(\vec{r},\vec{r'}) = \delta(\vec{r},\vec{r'}),
\end{equation}
那么根据上面的广义格林定理，就会得到如下形式
\begin{equation}
\int_V uLv d\tau = \int_V u \delta(\vec{r},\vec{r'}) d\tau
= u(\vec{r'}) = \int_V vLu d \tau = \int_V vf d\tau,
\end{equation}
这就是格林函数方法。

换一个角度，我们也可以这样看这个问题：方程原本是$Lu=f$，但我们通过构造格林函数$v$，成功地构造了算子$L$的逆算子，作用在$f$上，还原出了$u$，
\begin{equation}
u = \int_V vf d\tau.
\end{equation}
可以证明，格林函数$v(\vec{r},\vec{r'}) = v(\vec{r'}, \vec{r})$。

\paragraph{格林函数的数值方法}
求解两点边值问题，可以用格林函数方法，比如下面的常微分方程\cite{Hjorth-Jensen}
\begin{equation}
u''(x) = f(x), ~~ u(0)=u(1)=0.
\end{equation}
可以假设$u'(0)=1$(影响归一化因子)，然后用龙格库塔方法从左往右推，得到$u_{\rightarrow}(x)$。
也可以假设$u'(1)=1$，然后从右往左推，得到$u_{\leftarrow}(x)$。
那么我们可以定义格林函数
\begin{equation}
G(x_1, x_2) = u_{\rightarrow}(x_<) u_{\leftarrow}(x_>)
= \left\{
\begin{aligned}
& u_{\rightarrow}(x_1) u_{\leftarrow}(x_2), ~~ x_1<x_2 \\
& u_{\rightarrow}(x_2) u_{\leftarrow}(x_1), ~~ x_1 \geq x_2 \\
\end{aligned}
\right.
\end{equation}
显然$G(x_1,x_2)=G(x_2,x_1)$。
那么有
\begin{equation}
\frac{d^2}{d x^2_1}G(x_1,x_2) = \delta(x_1-x_2) [ u_{\rightarrow}(x_2) u'_{\leftarrow}(x_1) - u'_{\rightarrow}(x_1) u_{\leftarrow}(x_2) ].
\end{equation}
定义
\begin{equation}
W(x) = u_{\rightarrow}(x) u'_{\leftarrow}(x) - u'_{\rightarrow}(x) u_{\leftarrow}(x),
\end{equation}
则有
\begin{equation}
\frac{d^2}{d x^2_1}G(x_1,x_2) = \delta(x_1-x_2) W(x_2).
\end{equation}
根据广义格林定理，
\begin{equation}
\int^1_0 u(x) G''(x,x') dx = \int^1_0 G(x,x')f(x)dx,
\end{equation}
即
\begin{equation}
u(x')W(x') = \int^1_0 G(x,x')f(x)dx.
\end{equation}
所以有
\begin{equation}
u(x) = W(x)^{-1} \int^1_0 G(x,x')f(x')dx'.
\end{equation}
就这样，我们通过构造和计算格林函数（用到龙格库塔），把问题转化为了一个简单的数值积分。

\section{应用：三维球对称势场的薛定谔方程}
在球对称势场 $V(r)$ 中，定态薛定谔方程为
\begin{eqnarray}
( - \frac{ \hbar^2 }{ 2m} \nabla^2 + V(r) ) \psi ( \vec{r} ) = E \psi( \vec{r} ),
\end{eqnarray}
假设在球坐标下，体系本征态可以写作
\begin{eqnarray}
\psi( \vec{r} ) = R(r) Y^m_l ( \theta, \phi ),
\end{eqnarray}
则有
\begin{eqnarray}
- \frac{ \hbar^2 }{ 2m } ( \frac{1}{r^2} \frac{ d }{ dr } r^2 \frac{ d }{dr } - \frac{ l(l+1) }{ r^2} ) R(r) + V(r) R(r) = E R(r).
\end{eqnarray}
为了简化上式的形式，可以定义 $R(r) = u(r) / r $，则有
\begin{eqnarray}
- \frac{ \hbar^2 }{ 2m } \frac{ d^2 }{ dr^2 } u(r) + ( V(r) + \frac{ l(l+1) }{ r^2} \frac{ \hbar^2}{2m} ) u(r) = E u(r).
\end{eqnarray}
上式还可以进一步化简，定义
\begin{eqnarray}
\kappa &\equiv& \frac{ \sqrt{-2mE} }{ \hbar },  \\
\rho &=& \kappa r, 
\end{eqnarray}
则（43）式变为
\begin{eqnarray}
\frac{ d^2 u }{ d \rho^2 } = [ 1 - \frac{ V(r) }{ E } + \frac{ l(l+1)}{ \rho^2 } ] u(\rho).
\end{eqnarray}
有时这个形式更方便：
\begin{eqnarray}
\frac{ d^2 u }{ d \rho^2 } = [ 1 + \frac{ 2mV( \rho / \kappa) }{ \kappa^2 \hbar^2 } + \frac{ l(l+1)}{ \rho^2 } ] u(\rho).
\end{eqnarray}
解这个问题的算法如下：

1. 确定边界条件。
在无穷远处，有 $u(\infty) = 0$。
由于 $R(r)$ 在 $r=0$ 处为有限值，所以 $u(r) = r R(r) $ 在 $r=0$ 处有 $u(0) = 0$，而 $u'(0)$ 可以任取，只会影响 $u(r)$ 的归一性。

2. 选取一个足够大的值 $r_{max}$，将 $[0, r_{max}]$ 等分为 $n$ 份。
设定初始值 $u(r_{max}) = 1E-5, u'(r_{max}) = 0$，用 numerov 方法，从右往左推知各节点的函数值。
若未遇到极值点，则在 $r_{max}/10$ 处停止；若遇到极值点，则在极值点处停止。
设停止处为 $r_m$，从右往左推得该处函数值为 $u_{<} ( r_m )$。

3. 用 numerov 方法，从 $u(0)=0, u'(0) = 1$ 推得 $r_m$ 处函数值 $u_{>}(r_m)$。

4. 将 $r_m$ 右侧所有函数值乘以 $ u_< (r_m) / u_> (r_m) $，即使得函数连续。

5. 调节本征值 $\rho_0$，使得 $ u (r_m) - u( r_m - h ) = u( r_m + h) - u( r_m )$，即使得一阶导数连续。

\subsection{氢原子能级}
氢原子的势场为
\begin{eqnarray}
V(r) = - \frac{ e^2 }{ 4 \pi \epsilon_0 r },
\end{eqnarray}
其中 $e = 1.6 \times 10^{-19} C$。
将 $V(r)$ 代入（44），得到\cite{Griffiths2005Introduction},
\begin{eqnarray}
\frac{ d^2 u }{ d \rho^2 } = [ 1 - \frac{ \rho_0 }{ \rho } + \frac{ l(l+1)}{ \rho^2 } ] u,
\end{eqnarray}
其中
\begin{eqnarray}
\rho_0 &=& \frac{ me^2 }{ 2 \pi \epsilon_0 \hbar^2 \kappa}.
\end{eqnarray}

根据（44,45,50），能量与 $\rho_0$ 的换算关系为
\begin{eqnarray}
 E = - [ \frac{ m }{ 2 \hbar^2 } ( \frac{ e^2 }{ 4 \pi \epsilon_0 } )^2 ] \frac{ 4 }{ \rho^2_0 }.
\end{eqnarray}

而理论上氢原子的本征能量为
\begin{eqnarray}
 E = - [ \frac{ m }{ 2 \hbar^2 } ( \frac{ e^2 }{ 4 \pi \epsilon_0 } )^2 ] \frac{ 1 }{ n^2 },
\end{eqnarray}
其中 $n$ 为主量子数，轨道角动量量子数 $l = n-1, n-2, ..., 0$。

所以，$\rho_0$ 的理论值为 $\rho^{th}_0 = 2, 4, 6, 8, ...$。
表中展示了计算出来的 $\rho^{cal}_0$ 的值，与理论值相符。

\begin{table}
\caption{ 氢原子能谱对应的 $\rho_0$ 计算值 $\rho^{cal}_0$，$\rho_0$ 的定义请见文中。
在计算中，我们选取 $\rho_{max} = 50$，将 [0, 50] 分为 10000 个区间。
}
\centering
\begin{tabular}{ccccc}
\hline \hline
	$l=0$		&	$l=1$		&	$l=2$		& $l=3$		&	$l=4$		\\
\hline
	2.00018	&				&				&				&				\\
	4.00063	&	3.99989	&				&				&				\\
	6.0025	&	6.0004	&	5.99992	&				&				\\
	8.00451	&	8.00154	&	8.00154	&	8.00054	&				\\
	10.006	&	10.0014	&	10.0014	&	10.0014	&	10.0013	\\
	12.0077	&	12.0013	&	12.0013	&	12.0013	&	12.0013	\\
	14.0098	&	14.0012	&	14.0013	&	14.0013	&	14.0013	\\
	16.0121	&	16.0011	&	16.0012	&	16.0012	&	16.0012	\\
	18.0147	&	18.001		&	18.0012	&	18.0012	&	18.0012	\\
	20.0175	&	20.0009	&	20.0011	&	20.0011	&	20.0011	\\
\hline \hline
\end{tabular}
\end{table}

\subsection{ Woods-Saxon 势}
Woods-Saxon 势定义如下\cite{BABrown2011Lecture},
\begin{eqnarray}
V(r) = V_o (r) + V_{so} (r) \vec{l} \cdot \vec{s} + V_c (r),
\end{eqnarray}
其中，$V_o(r)$ 为中心势
\begin{eqnarray}
V_o(r) = \frac{V_o  }{ 1 + e^{(r-R_o)/a_o} },
\end{eqnarray}
其中 $R_o = r_o A^{1/3}$，$V_o, r_o, a_o$ 为可调参数，势场形状为费米分布函数（这可能是经验的）。

$V_{so}(r)$ 为自旋轨道耦合势
\begin{eqnarray}
V_{so}(r) = V_{so} \frac{ r^2_{so} }{r} \frac{ d f_{so} (r) }{ dr } = - \frac{ r^2_{so} V_{so} e^{(r-R_{so})/a_{so}} }{ a_{so} r ( 1 + e^{(r - R_{so})/a_{so}} )^2 },
\end{eqnarray}
其中 $R_{so} = r_{so} A^{1/3}$，$V_{so}, r_{so}, a_{so}$ 为可调参数，这个函数形式是从 Dirac 的理论推导出来的\cite{BABrown2011Lecture}。

$V_c(r)$ 为库仑势
\begin{eqnarray}
V_c(r) = \left\{
\begin{aligned}
& \frac{ Z e^2 }{ R_c } ( \frac{3}{2} - \frac{ r^2 }{ 2 R^2_c } ), ~~ r \leq R_c, \\
& \frac{ Z e^2 }{ r }, ~~ r \geq R_c.
\end{aligned}
\right.
\end{eqnarray}
其中，$Z$ 为质子数，$R_c = r_c A^{1/3}$，$r_c$ 为参数。

\paragraph{共振态与散射态}
\cite{Thijssen2011Computational}

\section{作业}
1. $x(t)$满足微分方程
\begin{equation}
\frac{dx}{dt} = -x,
\end{equation}
且有初始值$x(0)=1$，试用欧拉法求得$[0,5]$范围内的解$x(t)$，将数值解与解析解$e^{-x}$画在同一张图上。

2. 对于平抛运动，其运动方程为
\begin{eqnarray}
\frac{d}{dt} x &=& v, \\
\frac{d^2}{dt^2} y &=& -g,
\end{eqnarray}
不妨设$v=1,g=10$，初始时刻$x=0,y=0$，
\begin{itemize}
\item [1)] 试用龙格库塔方法，求取质点$t=10s$时的坐标。
\item [*2)] 稍微修改程序，求得$y=-h=-5$时的$x$值。
\end{itemize}